Excitable Patterns in Active Nematics 



L. Giomi, 1 ' 2 L. Mahadevan, 1 B. Chakraborty, 2 and M. F. Hagan 2 

1 School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA 
2 Martin A. Fisher School of Physics, Brandeis University, Waltham, MA 02454, USA 

(Dated: November 18, 2010) 

We analyze a model of mutually-propelled filaments suspended in a two-dimensional solvent. The 
system undergoes a mean-field isotropic-nematic transition for large enough filament concentrations 
and the nematic order parameter is allowed to vary in space and time. We show that the interplay 
between non-uniform nematic order, activity and flow results in spatially modulated relaxation 
oscillations, similar to those seen in excitable media. In this regime the dynamics consists of nearly 
stationary periods separated by "bursts" of activity in which the system is elastically distorted and 
solvent is pumped throughout. At even higher activity the dynamics becomes chaotic. 



Colonies of motile microorganisms, the cytoskeleton and 
its components, cells and tissues have much in common 
with soft condensed matter systems (i.e. liquid crystals, 
amphiphiles, colloids etc.), but also exhibit phenomena 
that do not appear in inanimate matter. These unique 
properties arise when the constituent particles are ac- 
tive: they consume and dissipate energy to fuel inter- 
nal changes that generally lead to motion. When active 
particles have elongated shapes, as seen in cytoskeletal 
filaments and some cells, they undergo orient ational or- 
dering at high concentration to form liquid crystalline 
phases. The theoretical and experimental study of ac- 
tive materials has disclosed a wealth of emergent behav- 
iors, such as the occurrence of giant density fluctuations 
PQ , the emergence of spontaneously flowing states (some- 
times referred to as flocking) [2] , unconventional rheolog- 
ical properties [3] and new spatiotemporal patterns not 
seen in passive complex fluids [4] . Recently Schaller et al 
[5] observed collective motion, large-scale fluctuations in 
density and the degree of polar order, and swirling mo- 
tions in a motility assay consisting of highly concentrated 
actin filaments propelled by immobilized molecular mo- 
tors in a planar geometry. However, the range of be- 
haviors realizable in active materials and the connection 
between material characteristics and system dynamics re- 
main incompletely understood. In this letter, we show 
that active systems exhibit behaviors similar to those 
of excitable systems, showing relaxation oscillations that 
couple activity to spontaneous pulsatile flow with quies- 
cent periods in between, similar to biological pumps. We 
furthermore show that systems which undergo nematic 
ordering can give rise to large-scale swirling motions re- 
sembling those observed in Ref. [5j even in the absence of 
polar order. These behaviors were not seen in previous 
investigations because they arise only when the degree of 
nematic order is allowed to fluctuate in space and time. 

Our model consists of mutually-propelled elongated 
particles suspended in a solvent confined to two dimen- 
sions, as seen for example in recent motor- filament assays 
[5]. The dynamical variables in such a system are par- 
ticle concentrations c, the solvent flow field v and the 
nematic tensor Qij = S(niTij — with S the nematic 



order parameter and n the director field, all of which are 
allowed to vary in space and time. We assume that the 
suspension of rod-like active particles have length i and 
mass M in a solvent of density p so ivent- The total density 
of the system p = Mc + p so ivent is conserved, so the fluid 
is assumed to be incompressible. The total number of 
active particles is also constant, thus the concentration c 
obeys the continuity equation: 

d t c=-V-(j p +3 a ), (1) 

where j p and j a are respectively the passive and active 
contributions to the current density. The passive cur- 
rent density has the standard form jf = cvi — DijdjC 
where = D^Sij + DiQij is the anisotropic diffusion 
tensor, while the active current can be constructed phe- 
nomenologically to be of the form jf = —aic 2 djQij or 
derived from microscopic models [6]. Here the factor c 2 
reflects the fact that activity arises from interactions be- 
tween pairs of rods while the constant a\ describes the 
level of activity and is proportional to the concentration 
of motors and the rate of adenosine-triphosphate (ATP) 
consumption. 

The flow velocity obeys an active form of the Navier- 
Stokes equation 

pd t Vi = rjAvi - dip + djnj , (2) 

with rj the viscosity, p the pressure, and the active stress 
tensor r^- given by: 

Tij = —XSHij + QikHkj — HikQkj + oi,2C 2 Qij (3) 

Here V • v = 0, and the first three terms in Eq. Q repre- 
sent the elastic stress due to the liquid crystalline nature 
of the system, with Hij — —SF/5Qij the molecular ten- 
sor defined from the two-dimensional Landau-De Gennes 
free energy: 

F/K = J dA[\(V-Q) 2 + \(c*-c) trQ 2 + |c(trQ 2 ) 2 ] 

where K is the splay and bending stiffness in the one- 
constant approximation. At equilibrium, above the crit- 
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FIG. 1: (Color online) The velocity filed (left) and the direc- 
tor field (right) superimposed to a density plot of the concen- 
tration and the nematic order parameter for a 2 = 0.4 (top) 
and «2 = 3 (bottom). The colors indicate regions of large 
(green) and small (red) density and large (blue) and small 
(brown) nematic order parameter. For moderate values of a 2 
the flow consists of two bands traveling in opposite directions 
with the director field is nearly uniform inside each band. For 
large ot<z the flow is characterized by large vortices that span 
lengths of the order of the system size and the director field 
is organized in grains. 



ical concentration c*, S = y^2tr Q 2 = y/i — c* /c con- 
sistent with hard-rod fluid models where the isotropic- 
nematic (IN) transition is driven only by the concentra- 
tion of the nematogens. The last term was first intro- 
duced in Ref. [8 and represents the tensile/contractile 
stress exerted by the active particles in the direction of 
the director field n with a 2 a second activity constant. 

Finally, the nematic order parameter Qij satisfies a 
hydrodynamic equation that can be obtained by con- 
structing all possible traceless-symmetric combinations 
of the relevant fields, namely the strain-rate tensor Uij = 
\{diVj + djVi), the vorticity tensor uj^ = \{diVj 
and the molecular tensor Hi« [7], so that 



djVi) 



[3t + v-V]Qi. 
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^Hij + XSuij + Qik^kj — ^ikQkj (4) 



where 7 is an orientational viscosity, and the additional 
terms on the right-hand side describe the coupling be- 
tween nematic order and flow in two dimensions, with 
A the flow- alignment parameter which dictates how the 
director field rotates in a shear flow and affects the flow 
and rheology of active systems [2j |3] . 

The dynamics of such an active nematic suspension 
is governed by the interplay between the active forcing, 
whose rate r" 1 is proportional to the activity parameters 
ql\ and 0^2, and the relaxation of the passive structures, 



the solvent and the nematic phase, in which energy is 
dissipated or stored. The response of the passive struc- 
tures, as described here, occurs at three different time 
scales: the relaxational time scale of the nematic degrees 
of freedom i? 2 /(7 _1 i\~), the diffusive time scale £ 2 /Do, 
the dissipation time scale of the solvent pL 2 /^, where 
L is the system size. While the presence of three dimen- 
sionless parameters makes for a very rich phenomenology, 
we temporarily assume that the three passive time scales 
are of the same magnitude r p . When r a ^> r p , the active 
forcing is irrelevant and the system is just a traditional 
passive suspension. On the other hand, when r a ~ r p , 
the passive structures can balance the active forcing lead- 
ing to a stationary regime in which active stresses are ac- 
commodated via both elastic distortion and flow. Finally, 
when r a <C r p the passive structures will fail to keep up, 
leading to a dynamical and possibly chaotic interplay be- 
tween activity, nematic order and flow. In the rest of the 
paper, we quantify these different regimes. We first make 
the system dimensionless by scaling all lengths using the 
rod length i, scaling time with the relaxation time of the 
director field r p = £ 2 / (7 -1 if ) and scaling stress using the 
elastic stress a = Ki~ 2 . 

We start with a linear stability analysis of the hy- 
drodynamic equations about the homogeneous solution 
letting (p = {c, Q xx , Q xy , 2uj xy } = <p + 5<p(t) with 
<£o = {co, *So/2, 0? 0}? an d consider solutions with peri- 
odic boundary conditions on a square domain of the form 
Sip(jc,t) = ^mW exp[^f(na: + m|/)]. 

With this choice the linearized hydrodynamic equations 
reduce to a set of coupled of linear ordinary differen- 
tial equations for the Fourier modes <£> nm : dt^p n m — 
A nm Lp nrn . The first mode to become unstable is the 
transverse excitation (n, m) = (0, 1) associated with the 
block-diagonal matrix: 
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The mechanism that triggers the instability of the ho- 
mogeneous state relies on the coupling between local ori- 
entations and flow embodied in the matrix 601, whose 
diagonalization yields a lower critical value of a 2 ' 



a = 



47r 2 [2r7 + S 2 (l-A) 2 ] 
c§L 2 S (l - A) 



(5) 



We see that the critical value of the activity required for 
instability falls with increasing system size and increases 
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FIG. 2: (Color online) The hydro dynamic fields c, S, and 
v x along the y axis for a 2 = 0.4 The yellow region indicates 
the band visible in the top panels of Fig. [I] 

with the viscosity of the solvent. In general, shear flow 
causes the director field to rotate for A ^ 1 , which gener- 
ates elastic stress. For small activity, the elastic stiffness 
dominates and suppresses flow, while above we ob- 
serve collective motion. 

To go beyond this simple analysis, we numerically inte- 
grated the hydrodynamic equations on a two-dimensional 
periodic domain with initial configurations of a homoge- 
neous system whose director field was aligned along the 
x axis subject to a small random perturbation in density 
and orientation, with ql\ = r] = c* = Do = D\ = 1, 

A = 0.1 and L = 10. As predicted by the linear sta- 
bility analysis, at low activity the system relaxes to a 
stationary homogeneous state with v x = v y = and 
S = yl — c*/c. Above the critical value a^, the system 
forms two bands flowing in opposite directions. The solu- 
tion is constant along the flow direction (see Fig. [l]top) 
while the direction of the streamlines (in this case along 
the x direction) is dictated by the initial conditions. As 
shown in Fig. [2j the optima in the flow velocity corre- 
spond to the maximal distortion of the director field n. 
Variations in concentration c and the nematic order pa- 
rameter S are of order 2% with a minimum in £ at the 
center of a flowing band due to the balance between dif- 
fusive and active currents. We note that these density 
modulations parallel to the nematic director are charac- 
teristic of systems with nematic order, while the density 
bands orthogonal to the direction of alignment observed 
in Ref. [5 are consistent with polar order [2 . 

Upon increasing the value of the activity parameter 
0^2, the spontaneously flowing state evolves into a pul- 
satile spatial relaxation oscillator. Fig. [3] (top left panel) 
shows a plot of the x and y component of the flow veloc- 
ity in the center of the box for = 1.5. In this regime 
the dynamics consists of a sequence of almost station- 
ary passive periods separated by active "bursts" in which 
the director switches abruptly between two orientations. 
During passive periods, c and S are nearly uniform, there 
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FIG. 3: (Color online) Dynamics of active "burst" for — 
1.5. The flow velocity at the point x = ?/ = L/3is shown as 
a function of time over the course of a director field rotation 
(top left) and the director field is shown for the three labeled 
time points. Between two consecutive bursts the system is ho- 
mogeneous and uniformly aligned. During a burst, nematic 
order is drastically reduced in the whole system and the di- 
rector undergoes a distortion with a consequent formation of 
two bands flowing in opposite directions. After a burst, a sta- 
tionary state is restored with the director field rotated of 90° 
with respect to its previous orientation. 

is virtually no flow and the director field is either paral- 
lel or perpendicular to the x direction. Eventually this 
configuration breaks down and the director field rotates 
by 90° (see Fig. [3|. The rotation of the director field is 
initially localized along lines, generating a band of flow 
similar to those in Fig. [I] (top). The flow terminates 
after the director field rotates and a uniform orientation 
is restored. The process then repeats. 

This rotation of the director field occurs through a tem- 
porary "melting" of the nematic phase. As shown in Fig. 
|3j during each passive period the nematic order param- 
eter is equal to its equilibrium value So = yl-c*/c, 
but drops to ~ |So during rotation. The reduction of 
order is system- wide, but, as shown in the middle in the 
bottom-left panel of Fig. [3J is most pronounced along the 
boundaries between bands. Without transient melting, 
the distortions of the director field required for a burst 
are unfavorable for any level of activity. 

To illustrate the origin of the relaxation oscillations 
it is useful to construct a simplified version of the hy- 
drodynamic equations by retaining the minimal features 
responsible for the oscillatory phenomenon: the coupling 
between active forcing and the fluid microstructure and 
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FIG. 4: (Color online) (Left) The average nematic order pa- 
rameter (S) and the total shear stress a xy are shown over 
several bursts for oli — 1.5. (right) The frequency of bursts is 
shown as a function of oli. 

the variable nematic order embedded in the Landau-De 
Gennes free energy. This can be achieved by approxi- 
mating Q 2 X as a constant and u xx ~ 0. Then, letting 
u = — u xy and Q = Q X y, and dropping the coupling be- 
tween Qij and c^, Eqs. Q and Q can be expressed in 
Fourier space as: 

Q = aQ-bQ 3 -u (6a) 
ii = k 2 (aQ — u) (6b) 

where A: is a wave number of an arbitrary spatial mode 
and a and b are constants. Eqs. Q has the form of the 
FitzHugh-Nagumo model for excitable dynamical sys- 
tems. For a < \{2a + k 2 ) the system rapidly relaxes 
to a state characterized by a finite strain-rate that bal- 
ances the active stress: u = aQ = otyj (a — a)/b. For 
a > \(2a + /c 2 ), this state becomes unstable and the 
trajectory converges to a limit cycle with a frequency 
v ^ k 2 a. As anticipated, when the active and passive 
time-scales are comparable the active forcing is accom- 
modated by the microstructure leading to a distortion 
of the director field and a steady flow. However, when 
the active forcing rate is increased, the microstructure 
dynamics lag, resulting in oscillatory dynamics. Clearly, 
the full system exhibits a much richer behavior than that 
captured by Eqs. (|6|, but the qualitative trends persist. 
For instance, the increase in the slope of v shown in Fig. 
[4] is not simply associated with the excitation of a second 
spatial mode of larger wave-number and is more likely 
due to the dynamics of the concentration field, which 
has been ignored in Eqs. ([6]). 

When the activity a 2 is further increased, the sequence 
of low activity periods and bursts becomes more complex 
and eventually chaotic. We emphasize that the oscilla- 
tory dynamics and the chaotic flows discussed below do 
not occur in our equations unless the magnitude of the 
nematic order parameter is allowed to fluctuate. Fig. [I] 
(bottom) shows a typical snapshot of the flow velocity 
and the director field superimposed to a density plot of 
the concentration and the nematic order parameter re- 
spectively. The flow is characterized by large vortices 
with positions related to "grains" in which the director 



field is uniformly oriented. The grain boundaries are of 
the order of the system size and are the fastest flow- 
ing regions in the system. Thus the dynamics in this 
regime is characterized by sets of grains of approxima- 
tively uniform orientation that swirl around each other 
and continuously merge and reform. 

In summary, we analyzed the hydrodynamics of ac- 
tive nematic suspensions in two dimensions. By allowing 
spatial and temporal fluctuations in the nematic order 
parameter we have observed a rich interplay between or- 
der, activity and flow. Significantly, we find that allowing 
fluctuations in the magnitude of the order parameter S 
qualitatively changes the flow behavior as compared to 
systems in which S is constrained to be uniform. Finally, 
we note that both the flip-flop dynamics (Fig. |3| and the 
swirling motion (Fig. [I] bottom) resemble behavior ob- 
served in the motility assay experiments of Schaller et al. 
[5]. Our analysis suggests that these classes of patterns 
can emerge even in the absence of polar order. 
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